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ABSTRACT 

We  provide  new  results  enabling  robust  interferometric  image  reconstruction 
in  the  presence  of  unknown  aperture  piston  variation  via  the  technique  of 
Redundant  Spacing  Calibration  (RSC).  The  RSC  technique  uses  redundant 
measurements  of  the  same  interferometric  baseline  with  different  pairs  of  aper¬ 
tures  to  reveal  the  piston  variation  among  these  pairs.  In  both  optical  and 
radio  interferometry,  the  presence  of  phase  wrapping  in  the  measurements  is 
a  fundamental  issue  that  needs  to  be  addressed  for  reliable  image  reconstruc¬ 
tion.  In  this  paper,  we  show  that  these  ambiguities  affect  recently-developed 
RSC  phasor-based  reconstruction  approaches  operating  on  the  complex  visi¬ 
bilities,  as  well  as  the  traditional  phase-based  approaches  operating  on  their 
logarithm.  We  also  derive  new  sufficient  conditions  for  an  interferometric  ar¬ 
ray  to  be  immune  to  these  ambiguities  in  two  different  senses:  immunity  up 
to  an  image  shift  in  the  reconstruction,  and  absolute  immunity.  We  show  the 
implications  of  these  results  for  imaging  via  phase  closures  and  extend  existing 
results  involving  the  classical  three-baseline  closures  to  generalized  closures. 
Furthermore  we  show  that  absolute  immunity  is  conferred  upon  arrays  whose 
interferometric  graph  satisfies  a  certain  loop-free  condition.  We  specify  this 
condition  and,  for  cases  in  which  this  condition  is  not  satisfied,  we  provide  a 
simple  algorithm  for  identifying  those  graph  cycles  which  prevent  its  satisfac¬ 
tion.  Finally  we  apply  this  algorithm  to  diagnose  and  correct  a  member  of  a 
pattern  family  popular  in  the  literature. 
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1  INTRODUCTION 

Optical  interferometry  is  a  multi-aperture  imaging  technique  which  is  attracting  increasing 
interest  in  the  astronomical  and  remote-sensing  communities.  The  appeal  of  this  technique 
is  primarily  due  to  the  high  resolution  it  affords  relative  to  single- aperture  imaging.  Namely, 
the  angular  resolution  of  a  single  aperture  is  limited  by  diffraction  to  where  A  is  the 
wavelength  of  the  light,  and  D  is  the  diameter  of  the  aperture.  On  the  other  hand,  the 
achievable  angular  resolution  of  an  interferometer  is  instead  given  by  A  ,  where  Bmax 
is  the  maximum  spatial  separation  of  any  two  telescopes  in  the  array.  Therefore  with  in¬ 
terferometry  one  can  achieve  the  same  high  resolution  offered  by  an  extremely  large  (and 
often  prohibitively-costly)  telescope  by  interfering  light  from  several  telescopes  of  practical 
size.  Optical  interferometers  image  a  scene  by  sampling  the  2D  Fourier  Transform  of  the 
scene.  Several  excellent  surveys  exist  which  describe  the  concept  of  interferometry,  including 
Labeyrie  et  al.  (2006),  Glindemann  (2011).  Each  pair  of  telescopes  measures  a  single  angular 
spatial  frequency  of  radians,  where  b  is  the  vector  difference  of  the  telescope  positions, 
which  is  known  as  a  baseline.  For  an  array  of  N  apertures,  the  data  set  then  consists  of  all 
(^)  such  measurements. 

A  principal  challenge  in  interferometry  is  variation  in  the  complex  gains  among  the  mul¬ 
tiple  apertures  of  the  interferometer.  In  radio  interferometry,  this  variation  can  arise  from 
differences  among  the  analog  components  of  the  antenna  elements  (e.g.  cable  length  differ¬ 
ences)  in  the  array.  In  optical  interferometry,  atmospheric  turbulence  distorts  the  wavefronts 
arriving  at  each  telescope  aperture  so  that  their  effective  path  lengths  from  the  target  (or  op¬ 
tical  pistons )  are  altered  by  a  random,  non-uniform  amount.  For  simplicity,  we  will  heretofore 
refer  to  such  aperture-specific  phase  variation  as  optical  path  difference  (OPD),  regardless 
of  its  source.  As  a  result  of  OPD,  the  Fourier  component  measured  by  apertures  i  and  j  is 
given  by  =  | yl3 \ e^°ij .  in  which  6i3  =  ( dij  -\-S4>ij)  mod  2n,  and  9tJ  is  the  true  Fourier  phase 
at  spatial  frequency  *yL,  and  Stfij  represents  the  OPD  observed  between  apertures  i  and  j. 
One  approach  to  eliminate  OPD  is  to  form  triple  products  of  the  Fourier  components  along 
the  sides  of  a  baseline  triangle  (e.g.  bi2,  b23,  and  b3i).  Note  that  the  OPD  cancels  in  these 
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products  and  hence,  like  the  Fourier  magnitudes,  these  so-called  bispectra  are  OPD-invariant 
observables.  However,  for  a  non-redundant  array  with  distinct  baselines,  recovery  of  the 
Fourier  phases  from  the  bispectra  phases  is  an  ill-posed  problem  since  there  are  only  (^.f1) 
independent  bispectra  (Readhead  et  al.  1988).  Successful  bispectra-based  image  reconstruc¬ 
tion  remains  feasible  in  spite  of  this  ill-posedness  (see  e.g.  Thiebaut  (2013),  Besnerais  et  al. 
(2008)),  but  in  this  case  prior  constraints  (e.g.  on  the  image  support)  must  be  enforced  to 
regularize  the  reconstruction. 

An  alternative  and  intrinsically  well-posed  approach  to  prior-regularized  reconstruction 
is  to  use  baseline  redundancy  to  explicitly  solve  for  OPD  variation;  an  array  with  baseline 
redundancy  contains  repeated  instances  of  the  same  baseline  involving  distinct  aperture 
pairs.  Since  Fourier  phases  can  be  assumed  to  be  equal  for  all  repeated  baselines,  an  observed 
difference  amongst  their  corresponding  measurements  exposes  the  contribution  of  the  OPD. 
This  idea  of  using  redundant  arrays  to  calibrate  out  OPD  variation,  known  as  redundant 
spacing  calibration  ( RSC ),  was  developed  in  works  such  as  those  by  Arnot  et  al.  (1985) 
and  Greenaway  (1990).  In  recent  years,  innovation  in  optical  technology  has  engendered  a 
revival  of  interest  in  the  RSC  technique.  The  simultaneous  (or  Fizeau- style)  measurement  of 
fringes  on  a  common  focal  plane  has  long  been  a  popular  method  of  acquiring  many  baseline 
measurements  in  an  economical  manner.  However,  the  Fizeau  method  had  been  incompatible 
with  RSC  techniques  since  the  fringes  formed  by  each  set  redundant  baselines  would  alias 
on  the  focal  plane.  An  elegant  solution  to  this  problem  was  proposed  by  Perrin  et  al.  (2006). 
This  work  developed  the  idea  of  segmenting  the  entrance  pupil  of  a  single  telescope  into 
an  RSC  arrangement  of  sub-pupils  from  which  the  light  was  then  coupled  via  single  mode 
fiber  to  a  non-redundant  exit  pupil,  thereby  permitting  unambiguous  and  simultaneous 
fringe  detection  for  an  RSC  array.  A  reconstruction  algorithm  for  this  architecture  was  then 
proposed  in  Lacour  et  al.  (2007).  Even  more  recently,  RSC  has  been  implemented  as  the 
calibration  scheme  of  choice  for  several  radio  interferometers:  the  Donald  C.  Backer  Precision 
Array  for  Probing  the  Epoch  of  Reionization  (PAPER)  in  South  Africa  (see  Ali  et  al.  (2015)), 
the  MIT  Epoch  of  Reionization  (MITeOR)  in  the  United  States  (see  Zheng  et  al.  (2014)), 
and  the  Ooty  Radio  Telescope  (ORT)  in  India  (see  Marthi  &  Chengalur  (2014)). 

As  will  be  shown  below,  N  —  3  independent  redundant  relations  are  required  for  unique 
determination  of  atmosphere  and  Fourier  phases  -  a  condition  we  will  refer  to  as  critical 
redundancy.  An  oft-cited  drawback  of  the  RSC  approach  is  that  it  reduces  the  number  of 
unique  spatial  frequencies  measured  by  the  interferometer.  However,  as  Figure  1  illustrates, 
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Figure  1.  Fraction  of  redundant  baselines  required  for  critical  redundancy  vs.  aperture  count 


the  fraction  of  distinct  uv-samples  sacrificed  for  critical  redundancy  becomes  increasingly 
negligible  as  the  number  of  apertures  in  the  array  increases.  Nevertheless  the  RSC  technique 
presents  other  challenges  which  must  be  overcome  for  reliable  imaging.  Central  among  these 
challenges  is  the  problem  of  integer  phase  ambiguities  which  arise  from  the  fact  that  the 
interferometric  phase  is  only  known  modulo  2n.  In  this  paper,  we  describe  these  ambiguities 
and  how  they  can  be  mitigated  using  a  combination  of  lattice  theory  algorithms  and  careful 
array  design.  We  will  see  that  these  ambiguities  have  a  fundamental  presence;  namely,  they 
exist  whether  the  calibration  strategy  works  with  complex  visibilities  (which  we  call  the 
phasor  approach)  directly,  or  their  respective  logarithms  (which  we  call  the  phase  approach). 

The  paper  is  organized  as  follows.  In  Section  II,  we  review  previous  work  on  the  integer 
ambiguity  problem,  and  discuss  its  presence  in  both  phasor  and  phase  approaches.  We  also 
provide  new  mathematical  conditions  for  an  aperture  pattern  to  be  wrap-invariant  in  the 
sense  that  it  is  immune  to  these  integer  ambiguities.  These  results  are  founded  upon  the 
well-known  Smith  Normal  Form  (SNF)  of  an  integer  matrix.  We  show  the  implications 
of  these  results  on  imaging  with  three  types  of  interferometric  observables:  the  baseline 
phase  measurements,  their  traditional  closure  phases,  and  generalized  closure  phases.  In 
Section  III,  we  relate  these  mathematical  conditions  to  conditions  on  the  aperture  pattern 
itself.  Namely  we  show  that  wrap-invariance  is  conferred  upon  arrays  satisfying  a  certain 
loop-free  condition.  As  an  illustrative  example,  we  diagnose  a  pattern  belonging  to  the 
popular  Y -pattern  class  and  remedy  it  to  be  loop-free.  Finally,  in  Section  IV,  we  show  that 
the  computationally-complex  SNF-based  approach  for  ambiguity  resolution  is  actually  not 
necessary  for  wrap-invariance  in  many  cases,  as  long  as  a  piston-dependent,  image  shift  can 
be  tolerated.  We  provide  physical  intuition  for  this  simpler  sufficient  condition  for  wrap- 
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invariance  in  terms  of  a  generalized  notion  of  the  well-known  concept  of  closure  phase. 
Finally  we  summarize  our  results  in  Section  V. 

2  PHASE  WRAPPING  AMBIGUITIES  IN  RSC  IMAGE 
RECONSTRUCTION 

2.1  The  Phase  Approach 

The  traditional  approach  to  RSC  reconstruction  operates  on  the  measured  baseline  phases 
(see  e.g.  Arnot  et  ah  (1985),  Greenaway  (1994)).  To  illustrate  the  approach,  let  us  consider  an 
interferometer  which  operates  at  a  wavelength  A  with  two  apertures  (say,  i  and  j )  separated 
by  a  vector  baseline  distance  of  hi3 .  In  the  absence  of  any  optical  path  difference,  the 
interference  pattern  formed  by  these  two  apertures  encodes  a  sample  of  the  object’s  Fourier 
Transform  at  spatial  frequency  'yh  Let  the  true  Fourier  phase  (which  we  will  refer  to  as 
object  phase),  measured  by  this  interference  pattern  be  denoted  as  0l3 .  The  measured  phase 
is  given  by: 


Pij  =  Oij  +  (j>j  -  fa  +  27re  (1) 

where  (pj  —  (pi  is  the  optical  path  difference  between  apertures  j  and  i,  and  e  is  unknown 
phase  wrap  integer  arising  from  the  fact  that  interferometric  phase  measurements  are  only 
known  modulo  2n. 

Consider  an  interferometric  array  which  simultaneously  makes  many  such  measurements 
amongst  its  N  apertures.  Suppose  that  of  the  array’s  (^)  baselines,  d  of  them  are  distinct. 
Further  suppose  we  have  a  solution  set  {(pi}  and  {0^}  for  these  equations.  Let  r,  denote  the 
vector  position  of  the  i-th  aperture.  As  noted  by  several  authors  (see,  e.g.  Wieringa  (1992)), 
we  can  obtain  another  valid  solution  set  simply  by  replacing  each  (pi  with  fa?  =  fa  +  epo  +  z-ri, 
and  each  9l3  with  9 =  0,3  —  z  •  (r,  —  r*),  for  arbitrary  <po  and  z.  Since  the  free  vector  z  is  a 
two-parameter  vector  representing  the  inherently-ambiguous  position  of  the  image  within  the 
Field-of-View  and  the  free  parameter  (p0  is  simply  a  scalar  piston  offset,  the  kernel  of  the  RSC 
system  is  three-dimensional.  Therefore  the  RSC  system  contains  d  unknown  distinct  object 
phases  and  N  unknown  aperture  pistons,  and  is  rank-deficient  by  at  least  3.  This  implies 
that  there  are  at  most  (d  +  N  —  3)  linearly-independent  equations  in  the  RSC  system,  and 
hence  at  most  N  —  3  redundant  relations  can  be  linearly-independent.  We  will  assume  for 
the  remainder  of  the  paper  that  our  array  contains  N  —  3  independent  relations.  Linder  this 
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Figure  2.  5- Aperture  Redundant  Array 


assumption,  we  can  then  solve  for  a  particular  solution  of  this  system  by  arbitrarily  setting 
two  object  phases  (whose  spatial  frequencies  are  not  co-linear)  and  one  piston  phase.  This 
particular  solution  will  then  differ  from  the  true  solution  by  a  phase  ramp  in  the  Fourier 
domain,  corresponding  to  an  image  shift  in  the  spatial  domain. 

As  an  example,  consider  the  simple  array  in  Figure  2.  There  are  (j)  =  10  baselines,  of 
which  4  are  redundant.  A  critical  array  of  5  apertures  would  have  2  redundancies.  Therefore 
this  array  possesses  more  redundancies  than  necessary  (call  it  strongly  redundant) ,  and  we 
anticipate  that  the  resulting  system  will  be  overdetermined. 

The  measurement  equations  associated  with  this  array  can  be  written  in  matrix  form: 


1  0  0  0  0  0  1 

0  1  0  0  0  0  0 

0  0  1  0  0  0  0 

0  0  0  1  0  0  0 

0  0  0  0  1  0  1 

0  0  0  1  0  0  0 

0  0-10001 
-10  0  0000 
0  0  0  0  0  1  0 

0  1  0  0  0  0  1 

Denoting  the  matrix  above  as  M, 


-10  0 
1  -1  0 

0  1  -1 

0  0  1 

0-10 
1  0  -1 

0  0-1 
0  1  0 

1  0  0 

0  0  0 


0 

0 

0 

-1 

0 

0 

0 

-1 

-1 

-1 


$ 12 
$23 
$34 
$45 
$13 
$25 
01 
02 
03 
04 
05 


=  0  +  27re 


we  can  write  this  system  in  compact  form  as: 


M 


$ 

0 


=  /3  +  2ne 


(2) 


(3) 


Since  the  interferometer  is  sensitive  only  to  optical  path  differences  amongst  its  apertures 
as  opposed  to  their  absolute  values,  one  of  the  0’s  (say,  0i)  can  be  set  to  zero  arbitrarily. 
Also,  as  noted  above,  we  can  set  two  of  the  object  phases  (say,  $i  and  $2)  arbitrarily.  After 
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removing  the  corresponding  three  columns  from  the  matrix  above,  we  are  left  with  a  full- 
rank  matrix  M.  If  the  wrap  vector  e  can  be  determined,  the  measurement  matrix  M  will 
admit  a  unique  solution,  and  the  RSC  phase  solution  is  in  hand.  We  hence  focus  our  analysis 
on  determination  of  the  wrap  vector. 

The  example  above  illustrates  that  the  general  phase  measurement  matrix  will  have  two 
sets  of  columns:  one  corresponding  to  the  object  phases,  and  one  corresponding  to  the  path 
differences.  Adopting  the  notation  of  Lannes  &  Anterrieu  (1999),  let  K  denote  the  subspace 
spanned  by  the  first  set  of  columns,  and  L  the  subspace  spanned  by  the  second  set. 

If  we  let  n  =  (i))  be  the  number  of  baselines  in  the  array,  the  phase  measurement  matrix 
M  will  be  of  size  n-by-(d  +  N  —  3).  For  a  strongly-redundant  array  like  the  one  in  the 
example  above,  the  column-space  K  +  L  of  the  matrix  will  clearly  not  span  Mn.  Therefore 
the  wrapped  measurement  vector  (3  will  not  in  general  fall  in  the  the  subspace  K  +  L  (and 
potentially  can  be  quite  far  from  it).  In  the  absence  of  measurement  noise,  we  can  unwrap 
these  measurements  by  identifying  those  integer  correction  vectors  e  for  which  y e  =  (3  +  2ne 
lies  in  K  +  L.  In  the  presence  of  noise,  on  the  other  hand,  the  unwrapped  vector  will  not 
generally  lie  in  K  +  L  (but  for  low-to-moderate  noise  will  be  in  the  vicinity).  Thus  we  search 
for  vector(s)  ye  which  are  as  close  to  K  +  L  as  possible  in  a  weighted  least-squares  sense 
(Lannes  &  Anterrieu  1999),  i.e.  we  search  for  the  vector 


Qrsc 

3 

tD 

1 

0 

j 

trsc  — 

=  argmine^f 

J 

4>rsc 

V 

0 

J 

(4) 

where  W  is  the  weighting  matrix1.  If  we  let  E  denote  the  phase  measurement  covariance 
matrix  and  set  W  =  S-1,  this  is  equivalent  to  searching  for  vectors  e  which  minimize  the 
projection  of  a  whitened  measurement  W^ye  =  W^(/3  +  2ne)  onto  the  space  ( K  +  L)fv  := 
/cer((W^M)T).  Specifically  we  seek  to  minimize: 


f(e)  =  \\PwW^f3  +  2ne)\\2  (5) 

where  is  a  matrix  representing  the  orthogonal  projection  from  M”  onto  ( K  +  L)fv. 
Letting  e'  =  — e,  we  can  rewrite  the  above  objective  function  as: 

/(e')  =  ||PVKW5(/0-2vre,)||2  =  |  |PwW2/3  -  2n¥>wW^e'\  f  (6) 

1  If  this  approach  successfully  unwraps  the  measurements,  the  least-squares  solution  associated  with  this  unwrapped  measure¬ 
ment  vector  will  be  the  Best  Linear  Unbiased  Estimator  (BLUE)  for  the  object  phase  vector  (see,  e.g.  Kay  (1993)) 
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Lannes  &  Anterrieu  (1999)  showed  that  this  optimization  problem  is  equivalent  to  the 
so-called  closest  vector  problem  in  the  theory  of  lattices.  We  will  define  a  lattice  L(Zn)  as 
the  set  of  points  generated  by  integer  combinations  of  the  column  vectors  of  a  matrix  L. 
Letting  P  =  P^W^,  our  optimization  problem  then  amounts  to  the  following:  Find  the 
lattice  point  in  P(Zn)  which  is  closest  to  P/3.  A  compact  representation  of  the  lattice  T  is 
given  by: 


{m^n—(d+N— 3)  j 

y  a HVi  I  Vcq  ez|  (7) 

where  {vj}  are  linearly-independent  and  together  form  a  basis  of  the  lattice.  Given  the  lattice 
basis,  several  algorithms  exist  for  finding  the  closest  lattice  point  to  a  specified  vector.  A 
popular  class  of  algorithms,  known  as  the  Sphere-Decoding  algorithms,  are  efficient  searches 
for  the  closest  lattice  point  within  a  hypersphere  of  a  certain  radius  centered  on  the  input 
vector  (see  e.g.  Agrell  et  al.  (2002)).  For  the  simulations  in  this  paper,  we  instead  use  the 
lower-complexity  Babai  Nearest  Plane  (Babai-NP)  algorithm  (Babai  1986).  For  lattice  bases 
which  are  nearly  orthogonal  (such  as  those  we  use  for  our  simulations),  this  algorithm  offers 
reliable,  albeit  not  guaranteed,  performance  in  practice. 

Suppose  we  have  found  a  basis  for  the  lattice  P(Zn),  and  we  have  solved  the  Closest 
Vector  Problem  for  a  given  measurement  vector  f3.  Let  b*  be  the  output  of  the  Babai 
Nearest  Plane  Algorithm  -  i.e.  it  is  the  lattice  point  which  is  the  closest  to  f3.  We  now  seek 
to  solve  for  the  wrap  vector  corresponding  to  this  lattice  point,  i.e.  we  seek  a  solution  to: 

b*  =  Pe  (8) 

Note  that  P  is  a  (weighted)  projection  matrix  and  thus  not  full-rank,  and  therefore  there 
will  be  inhnitely-many  solutions  to  this  equation.  The  Closest- Vector-Problem  algorithm  will 
provide  one  particular  solution  ep.  The  complete  set  of  solutions  is  then  given  by: 

e  Gp  T  G/j  (9) 

where  is  any  integer  vector  in  the  kernel  of  P.  Suppose  we  choose  one  such  vector  e ^ 
and  correct  our  phase  measurement  vector  accordingly.  The  corrected  phase  measurement 
vector  can  be  written  as: 

ft*  =  /3  +  27r(ep  +  e^)  (10) 
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Lemma  2.1:  G  K  +  L,  Ve^ 

Proof:  The  fact  that  G  ker( P^yW^)  implies  that  G  /cer(P^).  This  in  turn  implies 

that  W^e/,  G  im(W^M)  and  hence  that  G  im(M)  since  W?  is  invertible  by  construction. 

□ 

The  unwrapped  vector  f3*  therefore  differs  by  some  27re^  in  K  +  L  from  the  correctly- 
unwrapped  vector  y*,  i.e. 

/3*  =  y*  +  27re*  (11) 

We  now  examine  the  effect  of  this  residual  error  on  the  ultimate  least-squares  solution, 
which  is  easily  obtained  in  two  steps.  We  first  project  the  unwrapped  vector  onto  K  +  L. 
Noting  that  the  e*h  term  in  (3*  is  already  in  K  +  L,  we  obtain: 

ftc+L  =  W-i(I  -  P»-)wM*  =  y k+l  +  2lr<  (12) 

where  y*K+L  is  the  projection  of  y*  onto  K  +  L.  We  then  solve  the  system: 

M  \  =  4+i  (13) 

0 

where  we  have  chosen  to  work  directly  with  the  matrix  M  instead  of  M  for  purposes  of 
generality.  Since  M  is  rank-deficient  (by  3),  there  will  be  inhnitely-many  solutions  to  this 
system.  We  will  choose  a  solution  which  will  preserve  the  integrality  of  the  error  term  if 
possible,  thereby  yielding  a  final  RSC  solution  which  is  2ttc  away  from  the  true  solution  for 
some  integer  vector  c.  To  achieve  this,  we  rely  on  the  integer  matrix  decomposition  known 
as  the  Smith  Normal  Form,  which  is  described  in  the  following  Theorem: 

Theorem  (Smith  Normal  Form)  (Smith  1861):  Let  A  be  a  nonzero  m-by-n  integer  matrix 
with  rank  r.  There  exist  unimodular  (and  thus  invertible)  matrices  m-by-m  and  n-by-n 
matrices  U  and  V  respectively  such  that  the  matrix  product  D  =  UAV  is  a  diagonal 
matrix  whose  diagonal  entries  D,*  (the  so-called  elementary  divisors)  are  zero  for  i  >  r. 
Moreover,  the  product  of  the  elementary  divisors  is  the  greatest  common  divisor  (gcd)  of  all 
r-by-r  minors  of  A.  The  proof  of  this  Theorem  can  be  found  in  textbooks  such  as  Newman 
(1972).  □ 

Let  us  compute  the  Smith  Normal  Form  (SNF)  of  our  matrix  M: 

M  =  UmDmVm  (14) 
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where  the  r  diagonal  elements  {d{\  of  D M  are  the  elementary  divisors  of  M. 
We  can  now  re-write  Equation  (13)  above  as: 


D„V 


M  »  M 


1  Q* 


=  u 


M  t-’K+L 


We  can  choose  the  following  solution  to  Equation  (15): 


trsc  —  (dK+L 


(15) 


(16) 


where  Djy  denotes  the  pseudo-inverse  of  D.  The  resulting  error  is  then  clearly: 

eRSC  —  VM  D+Um  eh  (17) 

Lemma  2.2:  Let  u  =  Ujje*h  .  The  residual  wrap  error  e^sc  equals  0  mod  2tc  if  and  only  if 
mod  (ui,  di)  =  0 ,  Vi  ^  r.  The  proof  of  this  Lemma  is  an  adaptation  of  a  standard  proof  which 
can  be  found  in  most  textbooks  covering  linear  Diophantine  equations  (see,  e.g.  Newman 
(1972)).  □ 

From  this  Lemma,  the  following  Corollary  is  clear: 

Corollary  2.3  (Sufficient  condition  on  SNF  of  RSC  matrix  for  wrap -invariance ): 
If  the  elementary  divisors  of  the  measurement  matrix  M  corresponding  to  a  certain  aperture 
pattern  are  all  1,  the  RSC  solution  defined  by  trsc  is  immune  to  phase-wrapping  error. 

□ 


RSC  patterns  consisting  of  apertures  placed  randomly  on  a  Cartesian  grid  appear  to 
satisfy  this  sufficient  condition  with  high  probability.  We  conducted  a  simple  experiment  in 
which  15  apertures  were  randomly  placed  on  a  10-by-10  grid.  Out  of  256  placements,  66 
were  valid  RSC  patterns  (i.e.  possessed  at  least  critical  redundancy),  and  of  these,  only  2 
had  non-unity  elementary  divisors. 

It  is  noteworthy  that  in  cases  where  Corollary  2.3  holds,  we  may  be  able  to  compute  wrap- 
invariant  RSC  solutions  directly  (i.e.  without  actually  computing  the  SNF  of  M).  Namely, 
since  fi*K+L  is  in  the  subspace  K+L,  we  can  obtain  an  RSC  solution  by  solving  a  subset  of  the 
equations  in  Equation  (13).  Let  I  denote  the  indices  of  a  set  of  r  linearly- independent  rows 
of  M,  and  consider  the  sub-matrix  My  formed  by  selecting  only  these  rows  and  eliminating 
columns  correspond  to  two  arbitrarily-set  object  phases  and  one  arbitrarily-set  aperture 
phase.  Since  My  is  then  invertible,  we  can  form  a  possible  RSC  solution  as: 


TRSC  —  My  lfi*l,K+L 


(18) 
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The  resulting  error  6  in  the  final  RSC  solution  is  given  by: 

6  =  (19) 

The  error  will  consist  of  trivial  27T  errors  if  M71  contains  only  integer  elements.  It  is 
well-known  that  matrices  with  unitary  determinant  (so-called  unimodular  matrices )  have 
inverses  with  only  integer  elements.  We  thus  arrive  at  the  following  sufficient  condition  for 
the  solution  error  to  be  wrap-invariant  (i.e.  equal  to  zero  modulo  2n): 

Proposition  2.4  (Sufficient  condition  on  RSC  sub-matrix  for  wrap-invariance 
If  there  exists  a  unimodular  r-by-r  sub-matrix  of  M,  the  RSC  solution  TrSC  will  be  wrap- 
invariant. 

Note  that  no  such  unimodular  submatrix  will  exist  if  any  of  the  r  elementary  divisors 
of  M  are  not  equal  to  1,  since  in  this  case  the  gcd  of  all  r-by-r  minors  will  be  greater  than 
1  (see  Smith  Normal  Form  Theorem  above).  The  implications  of  this  Proposition  on  array 
design  will  be  made  clear  in  Section  3. 

The  SNF  has  been  been  applied  to  the  RSC  phase  problem  before  (Lannes  2003). 
Whereas  we  have  chosen  to  apply  SNF  directly  to  the  baseline  measurement  matrix,  the 
approach  taken  by  Lannes  (2003)  is  to  instead  treat  the  piston-invariant  phases  of  the  bis¬ 
pectra  (the  so-called  closure  phases )  as  the  fundamental  observables  from  which  the  object 
phases  can  be  inferred  via  the  relation: 

C  <,_*/  =  yd  +  2ned  (20) 

where  yd  are  the  wrapped  closure  phases,  ed  is  the  corresponding  wrap  vector,  and  C0^.c 
is  the  matrix  mapping  the  distinct  object  phases  in  the  array  to  closure  phases.  Lannes 
(2003)  hence  applies  the  SNF  to  the  closure  matrix  C0_>.c.  By  direct  analogy  to  Corollary 
2.3,  note  that  if  the  elementary  divisors  of  Coc  are  all  1,  then  the  pattern  is  wrap-invariant. 
Using  closure  phases  as  observables  can  be  advantageous  in  low-light  scenarios  in  which 
there  is  not  sufficient  SNR  in  a  single  atmospheric  coherence  time  to  reliably  measure  the 
baseline  phases.  To  overcome  this  low  per-frame  SNR,  atmosphere-invariant  observables 
such  as  the  bispectra  can  be  integrated  over  many  frames  to  build  sufficient  SNR,  and  their 
respective  closure  phases  used  as  reliable  phase  measurements.  Since  the  baseline  phases  are 
known  modulo  2i r,  the  linear  combinations  of  them  which  comprise  the  closure  phases  are 
also  known  modulo  27r.  In  order  to  relate  this  condition  to  Corollary  2.3,  let  us  first  define 
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Figure  3.  Distinction  between  spanning  tree  baselines  (thick,  solid)  and  loop  entry  baselines  (thin,  dotted) 


another  closure  matrix  Cm_>.c  which  instead  maps  the  phase  measurements  to  closure  phases. 
This  mapping  consists  of  equations  of  the  form: 


1/123  —  P 12  +  23  —  A 13  (21) 

where  7/123  is  the  closure  phase  associated  with  apertures  1,  2,  and  3,  and  the  are  the 
associated  baseline  phases  (see  Equation  (1)).  Of  the  Q)  possible  closure  phases,  at  most 
("k1)  can  be  linearly-independent  (see  e.g.  Readhead  et  ah  (1988)).  One  commonly-chosen 
set  of  such  linearly-independent  relations  consists  of  all  the  closure  triangles  involving  a  given 
aperture  A,  and  this  is  the  set  selected  by  Lannes  (2003).  Cm_>.c  is  therefore  an  ("“^-by-Q) 
matrix.  Lannes  (2003)  accordingly  provides  a  convenient  grouping  of  the  baselines  into  two 
categories:  (1)  spanning  tree  baselines  which  connect  aperture  A  to  all  other  apertures,  and 
(2)  loop  entry  baselines  which  provide  the  closure  for  these  spanning  tree  baselines.  This 
categorization  is  depicted  in  Figure  3. 

Given  this  categorization,  we  can  decompose  Cm_>c  into  corresponding  blocks  as: 


I 


(V)  J 


(22) 


where  Cm^c  contains  the  spanning  tree  contributions  to  the  matrix  (which  appear  in  mul¬ 
tiple  closures),  and  I^n-i\  is  the  ("9 1)_by-(n^1)  identity  matrix  representing  the  loop-entry 
contributions  (each  of  which  appears  in  only  one  closure).  The  following  property  follows 
from  this  block  form  expression: 

Lemma  2.5:  The  elementary  divisors  of  Cm_>.c  are  all  1. 

Proof:  Since  we  have  chosen  a  linearly-independent  subset  of  closure  relations,  r  =  rank( Cm4c 
(n^1) .  There  exists  a  r-by-r  minor  (namely,  I/n-i\)  which  is  equal  to  1.  Therefore  the  ged  of 
all  r-by-r  minors  is  1,  and  therefore  from  the  Smith  Normal  Form  Theorem,  all  elementary 
divisors  must  be  1.  □ 

Let  us  now  relate  C m_>c  to  the  matrix  C0_>.c  used  by  Lannes  (2003).  Recall  from  the 
discussion  of  bispectra  in  Section  1  that  the  closure  relations  eliminate  piston  differences 
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in  the  measurements  so  that  Cm_*.c  annihilates  the  subspace  L ,  i.e.  the  space  spanned  by 
the  columns  of  M  corresponding  to  <fi.  Defining  Mg  as  the  submatrix  of  M  containing  the 
columns  corresponding  to  6,  we  have 


9 

r  , 

ID 

Cm^cM 

l 

CUc  0 

1 

(23) 


where  CD_>c  =  Cm^.cMg.  C0_>c  is  an  (n“1)-by-(i  matrix  which  is  rank- deficient  by  two  2. 

Then  by  direct  analogy  to  Equation  (13),  we  can  obtain  valid  RSC  object  phase  solutions 
by  solving: 


0 

r  i 

6 

C,)WCM 

= 

CD^c  0 

=  y:,  +  2neld 

0 

0 

where  y*cl  is  the  true  unwrapped  closure  vector  and  2ne*h  d  is  the  residual  integer  wrapping 
error  vector  after  applying  the  Babai  NP  algorithm  to  solve  the  CVP  problem  associated 
with  matrix  Cm_*.cM0  and  2n ed.  Note  that  if  we  find  a  vector  6*  satisfying  Equation  (24),  it 
will  clearly  also  satisfy  the  relation  C o^,c0*  =  y*d-\-2/ne*h  d  of  Lannes  (2003).  Note  furthermore 
that  we  can  solve  the  equation  above  in  two  separate  integer-preserving  steps  of  the  form  of 
Equation  (16),  the  first  involving  the  SNF  decomposition  of  Cm_>c,  and  the  second  involving 
that  of  M.  Since  the  elementary  divisors  of  C.m_>.c  are  all  1  by  construction  (by  Lemma  2.5) 
and  hence  the  first  step  is  thus  integer-preserving,  wrap- invariance  again  amounts  to  whether 
or  not  all  elementary  divisors  of  M  are  1.  Therefore  we  have  the  following  Proposition  relating 
wrap  invariance  for  closure  measurements  to  that  for  raw  phase  measurements: 

Proposition  2.6  (Sufficient  condition  for  wrap -invariance  of  closure-based 
RSC):  If  the  elementary  divisors  of  the  phase  measurement  matrix  M  are  all  1,  then  the 
closure-based  RSC  solution  will  be  wrap- invariant.  □ 

We  remark  in  passing  that  although  the  preceding  analysis  was  presented  in  the  context 
of  the  traditional  three-aperture  closure,  it  applies  directly  to  the  case  of  closures  involving  an 
arbitrary  number  of  sides.  As  an  example,  consider  the  pattern  shown  in  Figure  4.  A  spanning 
tree  for  the  pattern  consisting  of  the  short  baselines  in  an  array  is  depicted.  Let  {<f>sp} 
denote  the  aperture  phase  differences  in  these  n  —  1  spanning  tree  baselines.  Note  that  all 


2  The  kernel  of  C 0— is  a  two-dimensional  subspace  of  the  three-dimensional  kernel  of  M.  To  see  this,  note  that  each  solution 
set  to  Equation  (20)  above  remains  valid  after  replacing  each  Oij  with  0 =  Oij  —  z  •  (r^  —  r^) 
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Figure  4.  Bootstrapping  phase  of  a  low-SNR  baseline  (green)  with  subset  (blue)  of  high-SNR  baselines  from  spanning  tree 
baselines  (black) 


aperture  phase  differences  in  the  array  can  be  expressed  as  linear  combinations  of  the  {0sp}. 
If  the  aperture  phase  differences  are  known  reliably  via  measurements  of  the  spanning  tree 
baselines,  we  can  use  these  measurements  to  cancel  the  aperture  phase  differences  in  all  other 
measurements  (of  which  one  example  is  shown  in  green).  The  idea  of  using  such  generalized 
closure  phases  (Martinache  2010)  is  indeed  the  mathematical  foundation  for  the  promising 
technique  known  as  baseline  bootstrapping  in  which  high-fidelity  phase  measurements  of 
several  high-SNR  baselines  (typically  the  short  baselines)  to  cancel  the  atmosphere  on  each 
lower-SNR  (and  hence  lower  fidelity)  baseline. 

Note  that  for  an  arbitrary  n-aperture  pattern,  there  will  in  general  be  n  —  1  spanning  tree 
baselines  and  thus  (")  —  (n  —  1)  =  (n“  )  generalized  closures,  each  involving  a  distinct  clos¬ 
ing  (or  loop-entry)  baseline.  Therefore  the  resulting  measurement  matrix  can  be  expressed 
exactly  as  in  Equation  (22)  above  and  hence  the  preceding  analysis  holds. 

While  this  section  has  considered  a  few  possibilities  for  phase  observables,  relating  mathe¬ 
matical  conditions  for  wrap  invariance  to  a  physical  condition  on  aperture  placement  is  more 
intuitive  when  considering  the  raw  phase  measurements  as  opposed  to  their  closures.  For 
this  reason,  for  the  remainder  of  the  paper  we  will  work  directly  with  the  baseline  phases  and 
their  associated  wrapping  errors.  In  particular,  we  will  begin  by  connecting  these  wrapping 
errors  with  analogous  ambiguities  in  recently-developed  phasor-based  approaches. 


2.2  The  Phasor  Approach 

Though  traditional  treatments  employ  the  phase  approach  of  the  previous  section  which 
operates  on  baseline  phases,  recent  papers  (e.g.  Marthi  &  Chengalur  (2014),  Liu  et  al.  (2010)) 
have  shown  that  a  non-linear  least  squares  (NLS)  approach  which  operates  at  the  phasor 
level  is  superior  in  accuracy.  Liu  et  al.  (2010)  developed  a  Gauss-Newton-type  NLS  solver 
and  showed  it  produced  unbiased  phase  estimates,  in  contrast  with  the  biased  ones  provided 
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by  the  phase  approach.  Marthi  &  Chengalur  (2014)  proposed  a  lower-complexity  alternative 
to  the  same  and  showed  that  it  achieved  performance  near  the  Cramer-Rao  Bound. 

The  Phasor  Approach  employs  the  following  measurement  model: 

Vij  =  gig*fij  +  nij  (25) 

where  Vt]  is  the  complex  visibility  observed  between  apertures  i  and  j,  gt  =  | gt  \  eJ^1  and 
dj  =  \dj\e^j  are  the  complex  gains  of  these  apertures,  ft]  is  the  true  complex  visibility 
measured  by  this  pair,  and  n  is  complex  measurement  noise.  Note  that  the  phase  difference 
between  gt  and  gj  is  simply  the  optical  path  difference  between  apertures  i  and  j  introduced 
in  the  previous  section.  Given  this  model,  NLS  approaches  attempt  to  find  a  set  of  complex 
phasors  {gi}  and  {/}  which  minimize  an  objective  function  of  the  form: 

A  =  EE'^ikh, -Mihm  -  stmtm  <26) 

i  j>i 

Minimization  of  A  with  respect  to  the  unknowns  (i.e.  distinct  object  and  antenna  com¬ 
plex  gains)  yields  the  following  conditions,  as  reported  by  Marthi  &  Chengalur  (2014)  in 
the  context  of  radio  interferometry,  and  by  Laeour  et  al.  (2007)  in  the  context  of  optical 
interferometry: 


9k  = 


h  = 


S j^k  Wkj  9jfkj  Vkj 

wkj\9j\2\fkj\2 
S j>k  9k9jVkj 


(27) 


(28) 


’Yhj>kWkj\9k\  \dj\ 

where  the  {fb}  are  the  true  complex  visibilities  of  the  distinct  object  phases  in  the  array. 

Due  to  the  circularity  of  these  definitions,  these  equations  must  be  solved  iteratively. 
Starting  from  an  initial  guess  for  all  phasors,  Equation  (27)  is  solved  to  obtain  a  better 
estimate  for  the  {gk}  and  then  these  {gk}  are  used  to  obtain  refined  estimates  of  the  {fb} 
through  Equation  (28).  In  the  next  iteration,  these  {fb}  are  used  to  further  refine  {gk},  and 
so  on. 

As  has  been  noted  before  (see  e.g.  Lannes  &  Anterrieu  (1999)),  there  are  strong  con¬ 
nections  between  phase-  and  phasor-based  approaches.  To  see  this,  let  z  be  the  vector  of 
products  {gj,g*f\i-j\}  which  minimize  A.  We  rewrite  Equation  (26)  as: 
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A  =  E5>«HWi  -  z‘i)(vi  -  4)H  (M) 

i  j>i 

Define  rl3  =  e^2nnijZij  for  an  arbitrary  integer  n^-  and  r  as  the  vector  containing  the  rl3 . 

Note  that  r  also  minimizes  A  since  the  rotations  {e?27rnu}  do  not  change  the  values  of  the 
residuals  in  A.  Hence  any  set  of  rotated  phasors  {&}  and  whose  products  produce 

the  vector  r  will  also  minimize  A.  Note  that  the  set  of  such  valid  phase  vectors  (i.e.  the 
concatenations  of  possible  {Ag*}  and  {A/|j_j|})  includes  the  set  of  phase  approach  solutions 
trsc  hi  Section  2  with  $*K+L  =  Ar  (where  Ar  is  the  vector  of  the  phases  of  the  complex  vector 
r).  In  other  words,  integer  ambiguities  present  in  the  phase  approach  do  not  disappear  in  the 
phasor  approach;  in  fact,  the  unwrapped  candidate  solutions  of  the  phase-based  approach 
correspond  to  local  minima  of  the  phasor-based  objective. 

Though  the  capacity  of  the  phasor  approach  to  produce  superior  accuracy  relative  to  the 
phase  approach  has  been  demonstrated,  the  former’s  convergence  issues  can  be  mitigated  via 
initialization  with  the  results  of  the  latter  (see  e.g.  Liu  et  ah  (2010),  and  Zheng  et  al.  (2014)). 
In  such  cases,  we  can  use  Corollary  2.3  from  the  previous  section  to  assess  the  viability  of 
the  phase-approach  solutions:  if  the  elementary  divisors  of  the  matrix  M  associated  with 
the  array  are  all  unity,  the  phasor-based  solutions  will  differ  from  the  true  solution  by  a 
multiple  of  2n.  Otherwise,  non-trivial  errors  are  possible  in  the  solution. 


3  IMPLICATIONS  OF  WRAP  AMBIGUITIES  ON  PATTERN  DESIGN 

In  this  section  we  use  the  mathematically-sufficient  conditions  for  wrap  invariance  from 
the  previous  section  to  show  that  aperture  patterns  whose  interferometric  graph  satisfies  a 
certain  loop-free  condition  are  wrap- invariant.  Here  we  define  the  interferometric  graph  in 
the  standard  way:  it  is  simply  the  graph  formed  by  connecting  the  array’s  apertures  (the 
nodes  of  the  graph)  with  edges  representing  the  array’s  baselines. 

This  condition  is  founded  on  the  sufficient  condition  in  Proposition  2.4  and  the  following 
definition  of  the  matrix  determinant.  This  definition  is  given  in  many  linear  algebra  texts 
(see  e.g.  Bretscher  (2001)). 

Definition  3.1  :  Suppose  we  have  an  n-by-n  matrix  A.  Define  a  pattern  as  a  selection 
of  n  entries  of  the  matrix  in  which  there  is  only  one  chosen  entry  in  each  row  and  one 
in  each  column  of  the  matrix.  Furthermore,  we  denote  a  pair  of  numbers  in  a  pattern  as 
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inverted  if  one  of  them  is  located  above  and  to  the  right  of  the  other.  Then  we  can  obtain 
the  determinant  of  A  by  summing  the  products  associated  with  all  patterns  with  an  even 
number  of  inversions  and  subtracting  the  products  associated  with  all  the  patterns  with  an 
odd  number  of  inversions. 

Our  goal  will  be  to  find  conditions  under  which  a  given  r-by-r  sub-matrix  Mi  contains 
only  one  pattern  with  a  non-zero  product,  in  which  case  the  determinant  will  be  ±1  from 
the  definition  above.  Consider  one  such  My  and  note  that  within  the  fully-connected  in¬ 
terferometric  graph  associated  with  M,  we  can  identify  a  sub-graph  G  containing  only  the 
measurements  in  My.  This  will  be  done  by  sequentially  identifying  those  matrix  entries 
which  must  be  part  of  a  pattern  with  a  non-zero  product.  Note  that  some  of  these  spe¬ 
cial  entries  from  the  matrix  My  can  be  identified  immediately.  Namely,  all  non- redundant 
measurements  contain  a  singleton  ±1  in  the  column  associated  with  their  object  phase.  All 
non-zero  patterns  must  clearly  contain  this  ±1  and  so  we  can  select  these  singleton  object- 
phase  entries  as  guaranteed  participants  in  a  non- zero  pattern.  Moreover,  all  measurements 
containing  a  leaf  node  (i.e.  a  node  with  a  single  connection)  in  G  contain  a  singleton  ±1  in 
the  column  associated  with  their  leaf  node.  All  non-zero  patterns  must  clearly  contain  this 
±1  as  well.  Thus  we  can  also  select  these  leaf  node  entries  as  guaranteed  participants  in  a 
non-zero  pattern. 

There  may  be  cascading  implications  of  such  singleton  measurements.  To  illustrate  this, 
consider  the  scenario  shown  in  Figure  5.  A  simple  RSC  array  is  shown  on  the  left.  A  subset 
of  the  baselines  in  one  possible  linearly- independent  sub-matrix  My  is  depicted.  Here  we 
intentionally  defer  selection  of  the  arbitrarily-set  phases  until  later  for  purposes  of  generality. 
A  simplified  depiction  of  the  matrix  My  is  shown  in  which  all  non-zero  entries  have  been 
colored  black  and  all  zero  entries  have  been  colored  white  for  simplicity.  In  Step  A  of  the 
reduction  process,  object  phase  9 5  is  selected  for  participation  (i.e.  its  matrix  entry  factored 
as  common  to  all  non-zero  patterns)  since  it  is  a  singleton  object  phase.  Its  corresponding 
row  (i.e.  row  5)  in  My  is  then  eliminated  from  participation  since  the  remaining  entries  in 
this  row  cannot  participate  in  a  pattern  (by  definition  of  a  pattern).  In  Step  B,  the  aperture 
6  entry  f>6  in  row  5  is  selected  since  it  has  become  a  leaf  node  in  the  pattern,  and  row  5  can 
then  be  eliminated.  Then  in  Step  C,  object  phase  64  is  then  selected  by  virtue  of  becoming  a 
singleton  object  phase,  and  row  4  is  then  eliminated.  This  selection/elimination  process  can 
be  repeated  beyond  the  steps  shown  in  the  Figure,  until  either  no  leaf  nodes  and  singleton 
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Figure  5.  Example:  Reducing  an  aperture  pattern  and  associated  matrix  to  identify  Persistent  Loop(s) 

object  phases  remain,  or  there  are  no  baselines  left  to  eliminate.  We  formalize  the  pattern 
reduction  process  in  Algorithm  1  below. 

Algorithm  1  Pattern  Reduction  Algorithm 

Require:  R  {where  R  is  the  set  of  the  baselines  corresponding  to  M/,  where  /  denotes  the 
indices  of  a  linearly-independent  subset  of  d  +  N  —  3  rows  of  M} 

while  \R\  >  0  do 

1.  Leaf  Nodes 

1.1  remove  any  remaining  baselines  containing  leaf  nodes  from  R 

1.2  add  the  associated  apertures  to  the  list  N 

2.  Singleton  Object  Phases 

2.1  remove  any  remaining  baselines  containing  singleton  object  phases  from  R 

2.2  add  the  associated  object  phases  to  the  list  O 

if  no  baselines  removed  in  the  current  iteration  then 
return  PERSISTENT 
end  if 
end  while 
return  LOOPFREE 

We  can  see  that  any  baseline  in  an  interferometric  graph  which  does  not  belong  to  a  loop 
will  be  eliminated  in  the  reduction  process,  and  its  corresponding  matrix  entries  factored. 
The  only  structures  in  the  graph  that  persist  after  this  reduction  are  sets  of  loops  with  a 
certain  property.  Namely  we  define  a  persistent  loop  set  as  a  set  of  loops  that  contains  at 
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least  two  instances  of  every  baseline  contained  in  the  set.  (A  set  can  consist  of  any  number  of 
loops,  including  one).  With  this  definition,  absolute  invariance  may  be  possible  if  the  graph 
of  the  redundant  baselines  does  not  contain  any  persistent  loop  sets.  Algorithm  1  returns 
PERSISTENT  if  persistent  loops  exist  and  LOOPFREE  if  the  pattern  is  completely  reduced 
and  therefore  free  of  persistent  loops. 

Note  that  in  the  latter  case,  we  will  have  eliminated  r  rows  from  My.  Since  each  baseline 
elimination  is  associated  with  object  phase  or  aperture  selection  from  distinct  columns,  and 
Mj  contains  r  +  3  columns,  there  will  be  exactly  three  extraneous  columns  not  involved 
in  the  reduction  process.  The  r-by-r  submatrix  obtained  by  selecting  the  non-extraneous 
columns  (i.e.  those  corresponding  to  the  selected  object  phases  and  leaf  nodes  in  O  and  N, 
respectively)  will  then  be  unimodular  by  virtue  of  having  a  single  non-zero  pattern  revealed 
by  the  reduction  process.  Having  ensured  the  existence  of  a  unit  r-by-r  minor,  we  have  hence 
confirmed  the  elementary  divisors  must  be  all  1,  and  hence  that  the  pattern  is  wrap- invariant 
(by  Corollary  2.3).  We  summarize  the  sufficient  condition  as  follows: 

Proposition  3.2  (Sufficient  conditions  on  aperture  pattern  for  wrap -invariance) 
Consider  the  graph  of  an  aperture  pattern  which  contains  d  distinct  baselines  and  any  set  of 
N —?>  linearly-independent  redundant  baselines.  If  this  graph  does  not  contain  non-persistent 
loop  sets  (in  the  sense  defined  above),  the  matrix  My  formed  by  these  independent  mea¬ 
surements  will  be  unimodular.  As  a  result  Proposition  2.4  from  Section  2  will  hold,  thereby 
guaranteeing  that  the  aperture  pattern  is  wrap- invariant. 

We  will  now  apply  Algorithm  1  to  the  example  pattern  shown  in  the  left  panel  of  Figure 
6.  This  pattern  belongs  to  one  of  the  more  popular  array  classes  in  the  interferometry 
literature:  the  so-called  Y-patterns  (see  e.g.  Arnot  et  al.  (1985),  Blanchard  et  al.  (1996), 
Labeyrie  et  ah  (2006),  Eastwood  et  ah  (2009),  Liu  et  ah  (2014)).  The  corresponding  spatial, 
or  UV,  sampling  is  provided  in  the  center  panel.  Algorithm  1  reduces  the  pattern  to  the 
persistent  loop  shown  in  the  right  panel. 

The  elementary  divisors  of  the  pattern’s  measurement  matrix  are  not  all  1;  they  are 
all  1  except  for  a  singleton  3  and  hence  detfWlj)  mod  3  =  0  for  all  choices  of  My.  To 
demonstrate  the  effect  of  the  resulting  phase  wrapping  on  reconstruction,  we  simulated 
noiseless  measurements  with  this  pattern  and  then  reconstructed  with  both  the  phase  and 
phasor  based  approaches.  The  results  are  shown  in  Figure  7.  The  upper  left  panel  shows 
the  true  image,  and  the  upper  right  panel  shows  the  projection  of  the  image  onto  the  space 
spanned  by  the  Fourier  basis  functions  measured  by  the  array  (i.e.  the  so-called  dirty  image). 
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Figure  6.  Example  of  a  pattern  containing  a  Persistent  Loop 

Y-Type  Aperture  Pattern  UV-Sampling  Persistent  Loop 

O 


Figure  7.  Reconstruction  Results  for  Y-Pattern  Example 
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The  lower  left  panel  show  the  reconstruction  result  with  the  SNF-based  phase  method,  and 
the  lower  right  panel  shows  the  same  for  the  phasor  method  using  the  updates  in  Equations 
(27)  and  (28).  Reconstruction  suffers  from  phase  wrapping  error  in  the  phase  case,  and  the 
corresponding  local-minimum  trap  in  the  phasor  case  as  discussed  in  Section  2.2.  The  closure 
phase  approach  yields  the  same  corruption  in  reconstruction,  as  the  elementary  divisors  of 
Coc  are  also  all  1  except  for  a  singleton  3. 

There  are  several  simple  ways  to  amend  this  pattern  so  that  it  is  wrap-invariant.  While 
the  most  intuitive  of  these  involve  moving  the  apertures  involved  in  the  persistent  loop  in 
Figure  6,  these  approaches  leave  gaps  in  the  UV-sampling  pattern.  An  alternate  approach 
that  preserves  the  UV-sampling  is  to  add  an  aperture  to  the  center  of  the  pattern  as  shown 
in  Figure  8.  This  results  in  additional  linear ly- independent  redundancies  colored  in  blue  and 
green,  respectively,  in  the  Figure.  These  additions  replace  baselines  in  the  persistent  loop, 
allowing  this  loop  to  be  broken.  With  wrap- invariance,  reconstruction  results  match  the  true 
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Figure  8.  Amended  Pattern 
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Figure  9.  Reconstruction  Results  for  Amended  Pattern 


Phase  Approach  (SNR  =  25  dB)  F 

3hasor  Approach  (SNR  =  25  dB) 

K 

image  in  both  the  phase  and  phasor  approaches  as  respectively  shown  in  Figure  9.  In  the  top 
row,  reconstruction  results  are  displayed  for  the  phase  (left)  and  phasor  (right)  approaches 
for  the  noiseless  case.  Analogous  results  for  an  SNR  of  25  dB  are  displayed  in  the  bottom 
row.  Here  we  define  SNR  as  the  ratio  of  the  phasor  magnitude  at  visibility  1  (i.e.  zero  spatial 
frequency)  to  the  standard  deviation  of  the  noise,  which  we  have  assumed  to  be  complex 
Gaussian  and  i.i.d.  across  spatial  frequency  for  this  simulation. 

4  CONDITIONS  FOR  WRAP-INVARIANCE  UP  TO  AN  IMAGE  SHIFT 

We  have  thus  far  seen  two  options  for  ensuring  wrap- invariant  RSC  solutions:  (1)  computing 
the  SNF  of  the  phase  measurement  matrix,  and  (2)  by  hireling  a  unimodular  submatrix  of  the 
measurement  matrix  M  via  selection  of  a  set  of  d  +  ( N  —  3)  linearly-independent  equations 
and  exclusion  of  two  object-phase  columns  and  one  aperture-phase  column.  Both  of  these 
options  may  become  computationally-burdensome  for  large  arrays,  as  in  those  under  current 
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consideration  with  Nap  fa  102  and  n  fa  104  baselines  (Zheng  et  al.  2014).  In  this  section  we 
develop  an  alternate  sufficient  condition  for  wrap  invariance  which  is  simpler  to  verify,  as 
it  does  not  require  the  computation  of  the  SNF  of  the  measurement  matrix.  This  condition 
applies  in  any  scenario  in  which  an  arbitrary,  piston-dependent  image  shift  can  be  tolerated. 

Note  that  the  spatial  frequencies  measured  by  an  array  can  be  represented  as  two-element 
vectors  of  the  form  (ujx,uy).  Let  X  be  the  d-by-2  matrix  containing  these  spatial  frequencies. 
Let  Mj  again  represent  the  submatrix  containing  a  linearly-independent  set  of  rows  of 
the  measurement  matrix.  Let  us  suppose  this  matrix  is  not  necessarily  unimodular  for  the 
scenario  in  question.  Let  B  =  M71.  Let  J  denote  the  indices  of  the  rows  of  B  which 
correspond  to  the  object  phases  associated  with  spatial  frequencies  in  X.  Finally,  let  e*h  be 
the  integer  error  vector  introduced  in  the  previous  section. 

Then  the  phase-wrap  error  will  manifest  itself  as  an  image  shift  if  and  only  if  this  error  is 
a  (modulo-27r)  phase  ramp,  i.e.  there  exists  a  2-element  shift  vector  z  and  an  integer  vector 
k  which  satisfy 

Bj(2ne*h)  -  2ttXz  =  2vrk  (30) 

Dividing  through  by  2ir  we  obtain  the  equation:  B  je  —  Xz  =  k.  Note  that  each  element 
of  B  j  can  be  expressed  as  some  rational  number  ^  where  q%  is  a  divisor  of  det(Mj).  Similarly 
we  first  assume  X  contains  rational  spatial  frequencies  with  greatest  common  denominator 
qx.  Then  we  can  multiply  through  by  the  least-common-multiple  (LCM)  of  (let (Mi)  and  qx 
to  obtain  a  system  of  equations  whose  coefficients  are  guaranteed  to  be  integer  (i.e.,  we  have 
a  linear  Diophantine  system).  Let  this  LCM  be  denoted  as  l.  Then  we  have,  after  rearranging 
terms,: 


/Xz  =  Z(Bje  -  k)  (31) 

We  now  wish  to  determine  conditions  under  which  there  exist  vectors  k  and  z  satisfying 
this  overdetermined  Diophantine  system.  Applying  the  Smith  Normal  Form  decomposition 
(see  Section  2.1)  to  the  matrix  /X  and  noting  that  rank(X.)  =  2,  we  have: 

Dx  =  Ux(/X)Vx  (32) 

where  Ux  and  Vx  are  unimodular  matrices  of  size  m-by-m  and  2-by-2,  respectively,  and 
Dx  is  a  rectangular  diagonal  matrix  whose  entries  are  zero  below  row  2. 
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If  we  left-multiply  Equation  (31)  by  U  on  both  sides,  we  obtain: 

/UxXz  = /Ux(Bje-k)  (33) 

Using  Equation  (32)  and  the  fact  that  V  is  a  unimodular  (and  hence  invertible)  matrix, 
we  can  then  write: 


DxV^z  =  /(UxBje  -  Uxk)  (34) 

We  then  arrive  at  the  following  proposition: 

Proposition  4.1:  A  pattern  is  invariant  to  integer  phase  ambiguities  up  to  an  image 
shift  if  the  matrix  UBj  contains  solely  integer  entries  below  row  2. 

Proof : 

Assume  the  matrix  UxBj  contains  solely  integer  entries  below  row  2.  We  re-arrange  the 
equation  above  so  that  it  reads: 

yDxV^z  —  UxB./e  =  —Uxk  (35) 

Let  v  =  jDxVx1z-UxBje.  Note  that  since  Dx  is  zero  below  row  2,  the  entries  of  v  below 
row  2  will  be  equal  to  those  of  (—UxBje),  which  are  integers  by  assumption.  Now  consider 
the  first  and  second  entries  of  v.  Let  f  be  the  vector  containing  the  fractional  parts  of  the 
first  two  elements  of  vector  UxBje,  and  let  A  be  the  invertible  matrix  consisting  of  the 
first  two  rows  of  yDxVx1.  Choose  z*  =  A_1f  so  that  the  fractional  part  f  is  annihilated, 
leaving  only  integer  elements  in  the  first  two  entries  of  v.  Hence  we  now  have: 


v  =  — Uxk  (36) 

with  v  ensured  to  contain  only  integer  elements.  Since  Ux  is  unimodular,  the  vector  k*  = 
— Uxxv  will  be  integral.  We  have  thus  found  a  pair  (z*,k*)  with  integer  k*  which  satisfies 
the  Equation  (34).  Since  Equation  (34)  is  related  to  Equation  (31)  via  a  unimodular  (and 
hence  invertible)  mapping  Ux,  invariance  is  hence  proven.  □ 

To  provide  an  intuitive  interpretation  of  this  condition,  we  revisit  the  notion  of  the 
generalized  closure  phase.  Let  us  dehne  a  generalized  closure  phase  as  an  integer  linear 
combination  of  object  phases  which  is  expressible  as  an  integer  linear  combination  of  phase 
measurements.  Note  that  in  order  for  a  combination  of  phase  measurements  to  satisfy  this 
condition,  the  corresponding  baselines  must  lie  on  a  closed  path;  otherwise  the  optical  path 
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terms  within  the  measurements  will  not  cancel  in  the  combination,  thereby  violating  our 
definition. 

Next  we  describe  the  function  of  the  operator  Ux  obtained  via  the  Smith  Normal  De¬ 
composition  (SNF).  The  SNF  matrices  Ux  and  Vx  associated  with  a  matrix  A  encode  the 
row  and  column  operations,  respectively,  necessary  to  convert  A  to  diagonal  form  Dx  In 
particular,  the  matrix  Ux  in  Equation  (34)  provides  the  integer  linear  combinations  required 
to  express  each  spatial  frequency  of  X  below  row  2  in  terms  of  the  spatial  frequency  vectors 
in  previous  rows  of  X.  That  is,  we  have: 

i 

X  u«xO.O  =  0  I  Vi  >  2  (37) 

3= 1 

where  X^.)  is  the  j'-th  row  of  X.  We  will  call  these  sums  the  canceling  combinations  of  the 
array’s  spatial  frequency  matrix  since  they  specify  how  to  cancel  each  new  spatial  frequency 
via  a  combination  of  previous  spatial  frequencies.  Now  recall  that  the  inverse  mapping  Bj 
maps  phase  measurements  f3  to  object  phases  9,  i.e.  we  have: 

9  =  Bj/3  (38) 

Left-multiplying  both  sides  by  Ux.  we  have: 

Ux0  =  UXB  jP  (39) 

The  integrality  of  UxBj  below  row  2  implies  that  integral  linear  combinations  (namely 
the  canceling  combinations)  of  object  phases  can  be  represented  as  an  integral  linear  com¬ 
binations  of  measurements.  In  other  words,  Proposition  4.1  is  equivalent  to  the  condition 
that  all  of  the  canceling  combinations  of  the  array’s  spatial  frequencies  correspond  to  a 
superposition  of  generalized  closure  phases  of  the  array. 

Note  that  verifying  the  sufficient  condition  in  Proposition  4.1  generally  requires  compu¬ 
tation  of  the  SNF.  However  this  SNF  computation  is  inherently  much  faster  than  the  SNF 
computation  in  Section  2,  since  here  we  are  computing  the  SNF  of  a  d-by-2  matrix  (relative 
to  a  n-by-(d  +  Nap)  in  the  previous  section).  Moreover,  assuming  the  unit  spatial  frequencies 
(0, 1)  are  (1,  0)  are  measured  by  the  array,  we  can  place  these  two  spatial  frequencies  at  the 
top  of  the  matrix  X  and  then  the  canceling  combinations  can  be  read  off  trivially  from  the 
elements  of  X. 

From  the  preceding,  we  can  obtain  a  couple  of  key  properties  of  arrays  that  satisfy 
Proposition  4.1,  which  we  will  call  the  image  shift  arrays  (ISA): 

(i)  Any  translation,  rotation,  or  scaling  of  an  ISA  will  also  be  ISA. 
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(ii)  Any  uniform  scaling  applied  to  either  of  the  dimensions  of  an  ISA  (i.e.  x  or  y)  will 
produce  an  array  which  is  also  ISA. 


Both  Properties  follow  from  the  fact  that  the  matrix  UxB  /  is  invariant  with  respect  to 
the  transformations  described.  Property  (ii)  is  of  interest  since  it  implies  that  the  ISA  prop¬ 
erty  is  preserved,  for  example,  when  an  array  is  altered  from  a  checkered  square  Cartesian 
grid  to  a  regular  hexagonal  grid  via  a  uniform  \/3-scahng  along  one  dimension. 

Finally,  it  is  interesting  to  note  that  it  is  highly-probable  that  a  randomized,  Cartesian- 
gridded  aperture  pattern  is  an  ISA.  To  demonstrate  this,  square  grids  of  various  sizes  were 
constructed  upon  which  apertures  were  placed  in  random,  uniform  fashion,  and  these  ran¬ 
dom  arrays  were  tested  for  the  ISA  property.  It  was  ensured  that  the  fraction  of  gridpoints 
occupied  by  an  aperture  was  constant  (at  ~  0.15)  for  all  grid  sizes.  This  constant  value  en¬ 
sured  that  random  placements  of  apertures  resulted  in  reasonably-sparse  arrays  that  were  at 
least  critically-redundant  with  high  probability.  The  grid  sizes  and  their  respective  aperture 
counts  are  shown  below.  We  plan  to  investigate  the  trend  shown  in  the  Figure  in  a  future 
paper. 


Grid  Size 

Area  (No.  of  gridpoints) 

No.  of  apertures 

12x12 

144 

22 

14x14 

196 

29 

16x16 

256 

38 

20x20 

400 

60 

24x24 

576 

86 

28x28 

784 

118 

The  results  of  the  experiment  are  shown  in  Figure  10. 
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Figure  10.  Empirical  ISA  Probability 


5  CONCLUSIONS 

We  have  examined  the  effect  of  phase  wrapping  in  Redundant  Spacing  Calibration  of  inter¬ 
ferometric  arrays.  In  particular  we  have  described  its  fundamental  presence  in  RSC  whether 
the  observables  considered  are  the  measured  baseline  phasors  or  their  phases.  Using  the 
Closest- Vector-Problem  formulation  due  to  Cannes  (2003),  we  have  developed  two  sets  of 
sufficient  conditions  for  an  array  to  be  immune  to  this  phase  wrapping.  From  this  analysis, 
we  can  outline  a  general  approach  for  ensuring  wrap  invariance  in  array  design.  If  an  arbi¬ 
trary  image  shift  in  reconstruction  can  be  tolerated,  it  suffices  to  verify  through  a  simple 
matrix  multiplication  that  canceling  combinations  of  the  array’s  spatial  frequencies  can  be 
mapped  onto  generalized  closure  phases  (c.f.  Section  4).  On  the  other  hand,  if  absolute  in¬ 
variance  is  required,  the  Smith  Normal  Form  of  the  measurement  matrix  can  be  computed. 
If  the  elementary  divisors  therein  arc  all  1,  the  pattern  is  wrap-invariant  and  the  Redundant 
Spacing  Calibration  solution  can  be  computed  as  shown  in  Section  3.  Furthermore,  we  show 
that  failure  to  meet  this  condition  amounts  to  the  existence  of  a  particular  kind  of  cycle  in 
the  interferometric  graph,  and  we  provide  an  algorithm  for  identifying  such  cycles  so  that 
they  can  be  removed. 
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